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By means of the variational approximation (VA) and systematic simulations, we study dynamics 
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^ and stability boundaries for solitons in a two-dimensional (2D) self-attracting Bose-Einstein con- 

's 

Q densate (BEC), trapped in an optical lattice (OL) whose amplitude is subjected to the periodic 

d time modulation (the modulation frequency, oj, may be in the range of several KHz). Regions of 
g 

I stability of the solitons against the collapse and decay are identified in the space of the model's 

C parameters. A noteworthy result is that the stability limit may reach the largest (100%) modu- 

1 lation depth, and the collapse threshold may exceed its classical value in the static lattice (which 
y—i corresponds to the norm of Townes soliton). Minimum norm Amin necessary for the stability of 
^ the solitons is identified too. It features a strong dependence on a; at a low frequencies, due to 
^ a resonant decay of the soliton. Predictions of the VA are reasonably close to results of the sim- 

ulations. In particular, the VA helps to understand salient resonant features in the shape of the 

o 

^ stability boundaries observed with the variation of u. 
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I. INTRODUCTION 



A challenging subject in the study of dynamical patterns in Bose-Einstein condensates 
(BECs) is the investigation of matter-wave solitons in multidimensional settings. Various 
routes leading to the creation of stable 2D and 3D solitons have been elaborated theoret- 
ically. The proposed approaches have much in common with their counterparts developed 
for the stabilization of 2D and 3D spatiotemporal solitons in nonlinear optics, see review [1] . 
However, thus far, only effectively one- dimensional (ID) matter- wave solitons have been cre- 
ated in experiments which used cigar-shaped traps to confine the condensate |2j (the actual 
shape of the soliton may be nearly three-dimensional, but it is intrinsically self-trapped only 
in the longitudinal direction). A phenomenon which impedes the stabilization of truly mul- 
tidimensional solitons is the occurrence of collapse in 2D and 3D condensates with attraction 
between atoms. As demonstrated in Refs. [31 HI [5], a universal method for the stabiliza- 
tion of matter-wave solitons in the multidimensional setting may be provided by optical 
lattices (OLs), i.e., periodic potentials which are induced, through interference patterns, by 
coherent laser beams illuminating the condensate in opposite directions (in principle, similar 
methods may be applied to a gas of polaritons, where an evidence of the BEC state was 
recently reported too [6]). The stabilization of 2D and 3D solitons is possible with the help 
of the fully-dimensional OL, whose dimension is equal to that of the entire space, and by 
low- dimensional lattices, whose dimension is smaller by one. The latter means quasi-lD and 
quasi-2D OLs in the 2D and 3D space, respectively [U [7j. OLs may also support solitons of 
a completely different type, namely, gap solitons in BEC (of any dimension) with repulsion 
between atoms. In this case, the solitons emerge due to the interplay between the repul- 
sive nonlinearity and the negative effective mass in parts of the linear bandgap spectrum 
generated by the OL [8]. Although gap solitons cannot represent the ground state of the 
respective system, they are dynamically stable objects [H]. The creation of effectively ID 
gap solitons, composed of ~ 250 atoms of ^^Rb, was reported in Ref. [TU], see also review 
[TT] . Multidimensional (chiefly, 2D) gap solitons [12j, and "semi-gap" solitons (which look 
like gap solitons in one direction, and regular solitons in the other [13]) were predicted too 
[12], but not yet observed in the experiment. Stable gap solitons, in ID and 2D settings 
alike, can be supported not only by periodic OLs, but also by quasi-periodic ones [H]. 

It is relevant to note that the presence of the lattice, defining a particular spatial scale 
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in the system (the OL period), introduces an effective nonlocahty. On the other hand, 
a nonlocahty can be directly induced by long-range interactions between atoms carrying 
dipolar moments. In this way, it has been demonstrated that the interplay between the local 
repulsion and long-range attraction may also support stable 2D solitons, both isotropic [15] 
and anisotropic [16j ones. 

Another theoretically elaborated stabilization technique relies on periodic time modu- 
lation of the effective nonlinearity coefficient, which may be provided by the Feshbach- 
resonance effect in a low-frequency ac magnetic field applied to the condensate. It has been 
demonstrated that this method can stabilize 2D fundamental solitons, but not 3D ones [T7j 
(the stabilization of 3D solitons and their bound complexes is possible under a combined 
action of this technique and a quasi-lD OL [20]; the application of the Feshbach-resonance 
technique to matter waves trapped in usual OLs [21], as well as in parabolic traps [H], was 
studied too). This approach to the creation of stable 2D matter- wave solitons followed the 
earlier elaborated mechanism providing for the stabilization of (2+l)-dimensional optical 
spatial solitons in the model of a bulk medium built as a periodic concatenation of layers 
with self- focusing and self-defocusing Kerr nonlinearity [18j. It should be noted, however, 
that very accurate simulations of the so stabilized 2D solitons in the framework of the respec- 
tive Gross-Pitaevskii equation (GPE) demonstrate that they may be subject to an extremely 
slow decay [23] . 

The methods using the periodic modulation of the nonlinearity belong to a broad class of 
techniques relying on the periodic management of solitons [23] ■ In terms of the BEG in one 
dimension, another relevant example representing this class is the lattice management in the 
model where the strength of the OL is subjected to periodic time modulation [251 ES] (in 
the experiment, this was used as a method to excite the condensate [27]). The main result 
of the analysis of the model including the lattice management in combination with the self- 
repulsive nonlinearity, reported in Ref. [2S], was identification of stability limits for several 
families of gap solitons in this setting (fundamental solitons, double-hump bound states, 
etc.). Also belonging to this general class are settings based on the periodic modulation 
in time of the strength of a parabolic potential trap which confines the condensate [28] , 
and periodic modulation of the transverse trapping frequency which supports cigar-shaped 
BEG configurations. In the latter case, the "management" may induce Faraday waves in 
the condensate, as was predicted theoretically [22] and demonstrated experimentally ^SU\ . 
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It is also relevant to mention that another variety of the "management" was proposed very 
recently, with the aim to stabilize 2D matter-wave solitons: driven Rabi oscillations in a 
mixture of two mutually interconvertible BEC species (different states of the same atom), 
with self-attractive and self-repulsive intra-species interactions [3T\ . 

II. THE MODEL AND OUTLINE OF THE PAPER 

A natural extension of the studies outlined above is the application of the lattice- 
management technique to multidimensional solitons. This topic includes several distinct 
problems, depending on the sign of the intrinsic nonlinearity. In the ID setting, the lattice 
management for solitons in the attraction model is not a fundamentally important issue, 
because stable ID solitons in the self-attractive condensate exist without any lattice [2]. 
However, the model with the intrinsic attraction poses a new problem in the 2D case, as, 
in the absence of the lattice, 2D solitons are unstable against the collapse, and, as said 
above, the simplest way to suppress the instability is provided by the use of the quasi- ID 
OL Therefore, it may be interesting to identify stability limits for 2D solitons in the 
attractive model with the quasi-lD lattice whose strength is subjected to the periodic time 
modulation; in particular, one may be wondering if the solitons may survive in the case when 
the modulation depth attains its maximum, so that the lattice periodically vanishes. This 
problem was a subject of recent work [32j, which was also dealing with collisions between 
the 2D solitons, that may move freely in the unconfined direction. It was found that the 
actual stability region is quite limited: typically, the solitons loose their stability before the 
modulation depth attains 50%. 

The objective of the present work is to construct soliton families and identify their stabil- 
ity boundaries in the 2D attractive model with the full 2D lattice subject to the time-periodic 
modulation. In the normalized form, the respective two-dimensional GPE for mean-field 
wave function \1/ [x, y, t) is 



where t is time, (x, y) are coordinates in the 2D plane (scaled so as to fix the OL period equal 
to vr), and Vq is the strength of the lattice, while e and u are the amplitude and frequency 
of its temporal modulation. Coefficient —1 in front of the nonlinear term in Eq. ([T]) implies 
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that the nonhnearity is attractive. In the case of the quasi-lD lattice [11[32], term cos(2?/) is 
dropped from the potential in Eq. ([T]) . Actually, Vq is measured in units of the recoil energy 
of atoms trapped in the OL. For a typical case of atoms of ^Li loaded in the OL with the 
period on the order of /im, a characteristic value of the scaled frequency, ~ 2 (see below), 
corresponds, in physical units, to the modulation rate on the order of several KHz. 

An experimental implementation of the model may be provided by periodic attenua- 
tion of the intensity of the laser beams that create the OL. Since this management mode 
cannot change the sign of the OL's amplitude, the modulation coefficient in Eq. ([T| is 
restricted to interval < £ < 2. We, chiefly, follow this restriction below. On the other 
hand, Eq. ([T]) with the modulated potential that formally corresponds to e oo, viz., 
—Vq cos(co't) [cos (2x) + cos (2?/)], may be realized as a superposition of several moving OLs 
[33] (a ID counterpart of that setting was introduced in Ref. [S])- In works [M] and [33] 
it was demonstrated that such potentials, if combined with the self-repulsive nonhnearity, 
give rise to anomalously narrow subdiffractive gap solitons, for which the usual second-order 
operator of the kinetic energy is, effectively, replaced by a fourth-order operator. On the 
other hand, a piecewise-constant time-periodic modulation mode was considered, for both 
periodic and quasi-periodic 2D potentials, in recent work [35j (also in the combination with 
the self-repulsive nonhnearity). The result was generation of a multi-soliton lattice, starting 
from an initial localized state. 

The rest of the paper is organized as follows. In Section III, we elaborate the variational 
approximation (VA) for the analysis of the dynamics of fundamental 2D solitons in the 
model based on Eq. ([T|. In Section IV, results of systematic direct simulations of the 
soliton evolution in Eq. ([T]) are reported, and compared with predictions of the VA. The 
results are summarized in the form of stability regions for the solitons in the model with the 
time-modulated OL. The most noteworthy flnding reported in Section IV is that, depending 
on uj, the solitons may remain stable up to the maximum (100%) modulation depth, in 
contrast to the shallow modulation that limited the stability of the solitons supported by 
the quasi-lD OL [32]. It is also worthy to mention a multi-resonance shape featured by the 
stability boundary in the {uj, e) plane, and an increase of the collapse threshold (deflned in 
terms of the soliton's norm) provided by the time- modulated OL, in comparison with the 
static lattice, i.e., the threshold may exceed the norm of the Townes soliton. The minimum 
norm necessary for the stability of the solitons, N^i-^, is identifled too, as a function of u 
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and e (this characteristic was considered earher only in the static model [36]). Due to the 
possibility of a resonantly enhanced decay of the soliton, A^min features a strong dependence 
on uj at low frequencies. The paper is concluded by Section V. 



III. THE VARIATIONAL APPROXIMATION 



A. Effective Lagrangian 

Variational methods have been quite useful in many problems of nonlinear optics and 
BEC [21 H HZl ESI EZl EH]. To apply the VA to the present model, we notice that Eq. g 
can be derived from Lagrangian L = Cdx, with density 



l + -cos(cjt) [cos(2x) + cos(2?/)] 1^1% (2) 



where the asterisk stands for the complex conjugation. Following Refs. [3] and [1], we adopt 
the isotropic ansatz for the soliton, 

M/ans (x, y, t) = Ait) exp (^icl){t) + '-bity - , (3) 

where = + y"^, and all variables A(t), (f>(t), b(t), and W(t) (amplitude, phase, radial 
chirp, and radial width, respectively) are real. Tractable, although rather cumbersome, 
variational equations may also be generated by an anisotropic generalization of ansatz (|3]), 
with different widths and chirps along axis x and y. However, subsequent analysis of those 
equations does not predict stable anisotropic solitons. 

The substitution of the ansatz in Eq. ^ and calculation of the integrals yield the effective 
Lagrangian, 



dt 2W^ AirW^ 



e , ■ 
1 + - cos(ci;t) 



2dt 2 ' ^ ' 

where N = ttA'^W'^. The first Euler-Lagrange equation following from effective Lagrangian 
(j4|, 5 (J Lcsdt) /5(f) = Q {5/5(j) stands for the variational derivative of the action functional), 
is tantamount to the conservation of the norm of the wave function, which is the single 
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dynamical invariant of Eq. ([T]). Indeed, the norm of ansatz ^ is 

\^.ns{x,y)fdxdy = TiA^W' = N. (5) 




The second Euler-Lagrange equation, 6 (/ Lf^^dt) /6b = 0, reduces to the well-known expres- 
sion for the chirp in terms of the time derivative of the width [[33 EH], b = {dW/dt). 
Using this relation, the next variational equation, which accounts to dLf^g/d {W"^) = [since 
Lagrangian Q does not contain dW/ dt] can be cast in the following final form, 

£W 1 - N/N, 



-Wo 



l + |cosM) Wexp{-W^), (6) 



df^ 

where A^max = Svr is the well-known VA prediction [2T| for the critical (maximum) norm in 
the 2D space, which separates collapsing solutions at > A^max, i-e., ones with W{t) —>■ 
0, A(t) ^ oo at t — > tcoUapse [for Vo = and initial conditions W(t = 0) = Wq and dW/dt(t = 
0) = 0, the collapse time predicted by Eq. (6) is tcoUapse = Wq (A^/A^max — l) ], and 
noncoUapsing ones at A^ < Nraax- 



B. Analysis of variational equation. 

The actual maximum value of A^, found numerically from Eq. ([T]) with Vq = (it gives 
the norm of the Townes soliton [39]), is 

A^ma. = 5.85 ^ 0.93iV^ax, (7) 

which characterizes the accuracy of the VA. For the condensate of ^Li atoms, a typical 
collapse threshold corresponds to the number of atoms ^10^. 

Equation ([6]) helps one to understand what may happen to the 2D soliton under the 
action of the weak "management", with e/2 <^ 1. First, for e = 0, Eq. ^ predicts a stable 
equilibrium position, which is given by a smaller root of equation 

<exp(-iyo^) = ^~^/^- (8) 

(the larger root gives an unstable solution). Then, the linearization of Eq. ^ (still with 
e = 0) yields the eigenfrequency of small oscillations around W = Wq, 



-0 - ^ . (9) 







For instance, in the example considered in the next section (see Fig. [5] below), with 1 — 
N/Nm^^ = 0.155 and Vq = 0.25, the relevant root of Eq. g is W^o ~ 0.71, and Eq. g yields 
ujq ^ 1.35. Keeping quadratic and cubic terms in the expansion of Eq. ^ in powers of 
w(t) = W(t)/Wo — 1 around W leads to a standard equation of driven nonlinear oscillations. 
In particular, in the near-critical situation, i.e., for 1 — A^/A^max ^ 1, this equation takes the 
form of 

— + IQVqw - 2AVow^ + AOVqw^ = -2eVo cos(cjt) - 2eVo cos (iut) w. (10) 

It predicts the lowest-order direct resonance when u is close to ujq, the parametric resonance 
at uj close to 2ujo, and higher-order resonances at = nujo, with n = 2,3,4,.. |10]. The 
resonances may help to stabilize the 2D soliton against the collapse, as, increasing the am- 
plitude of its intrinsic oscillations, the soliton spends less time in the "dangerous zone" with 
small width, that might be a starting point for the collapse. On the other hand, effectively 
stretching the soliton, the resonances may destabilize it against decay into radiation, at 
smaller values of A^. Both trends are observed in numerical simulations, as shown below. 
Comparison of predictions following from a numerical solution of full variational equation 
^ and direct simulations of Eq. ([T| is presented in the next section. 



IV. NUMERICAL RESULTS 
A. The simulation mode 

Systematic simulations of Eq. ([T]) were performed by means of the split-step method 
|42j , in the {x, y) domain of size 256 x 256 or 512 x 512. The simulations were run with the 
Gaussian initial configuration, 

^(x, y, 0) = Ao exp {-q [{x - xof + {y - y^f] ) , g > 0, (11) 

whose norm is = 'nA^/q [cf. ansatz Taking A^ < A^max, as well as A^ slightly exceeding 
Nraax [recall A^max, the critical value of the norm for the onset of the free-space collapse, is 
given by Eq. ([T])], it was quite easy to find stable solitons which keep their shape despite 
the temporal modulations imposed by the "management" , see generic examples in Figs. [T] 
and |2} The former figure presents comparison of the initial and final shapes of the soliton. 
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FIG. 1: (Color online) A typical example of a stable 2D soliton, as obtained from the numerical 
solution of Eq.Q with Vq = 0.65, e = 0.5, u> = 1.35, and initial configuration (11) with Aq = 1.39, 



q = 0.5, whose norm, N = 6.07, exceeds the collapse threshold in the static model, N^nax ~ 5.85. 
The panels display the density distribution, (x , y , t)\'^ , at t = (a) and t = 150 (b). 

while the latter one displays the evolution of the soliton's amplitude in the course of its 
self-adjustment to the stable shape. 

The soliton shapes displayed in Figs. [T] and [2] are confined, essentially, to a single cell of 
the OL (similar to results reported in previous works [31 El E])- Roughly the same shapes 
would be observed in a parabolic trapping potential; however, the difference is that, if the 
nonlinearity is too weak, the OL cellular potential cannot suppress the tunnel decay of the 
localized pulses [36J, therefore the decay is observed at smaller values of A^, as shown below. 

As said above, in the 2D equation with Vq = localized configurations with > A^max = 
5.85 suffer collapse, while those with A^ < A^max decay. The static OL stabilizes 2D solitons 
in the latter case, but it cannot arrest the collapse of initial localized states with A^ > A^max 
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FIG. 2: (Color online) (a) The evolution of the amplitude of the stable soliton, \A\ = 
1^' (x = xo,y = yo)\, for the same parameters as in Fig. [T] except for Vq = 0.5 (|^| has the same 
meaning in other figures). Curve f{t) shows the modulation function in Eq.Q, 1 + (e/2) sin (wt); 
(b) The shape of the soliton at t = 150. 

[3l m |5]. The numerical analysis of the model with the quasi-lD lattice potential subjected 
to the periodic time modulation did not reveal stable solitons with > A^max either [32]. As 
shown in Fig. [3| in the present model the periodic time modulation of the two-dimensional 
OL potential makes it possible to stabilize the solitons both at A^ < A^max and in some 
interval above A^max- The actual increase of critical norm is not large, but the very fact that 
the constraint A^ < A^max can be broken is an interesting result, as it has never been reported 
before, to the best of our knowledge. 
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FIG. 3: (Color online) The evolution of the soliton's amplitude, as obtained from the numerical 
solution of Eq. ([T| for e = 0.5 and u> = 1 and different values of norm N of initial configuration 



(11) (higher curves correspond to larger N). (a) The model without the optical lattice, Vq = 0. In 
this case, all configurations with N < A^max ~ 5.85 decay, while the ones with > A'^max suffer the 
collapse, (b) In the presence of the time-modulated optical lattice with Vq = 0.5, the 2D solitons 
may be stable, including some values of above A^max- In particular, the soliton is stable at 
AT = 6.01 « l.OSA^max, while it collapses for A^ = 6.09, at t = 1.49. If A^ is too small, A^ < A^min, 
the soliton gradually decays into radiation, see Fig. [Tjbelow. In panel (b), the curve for A^ = 3.14, 
which is slightly smaller than the corresponding value of A^miiD also shows the slow decay of the 
soliton. 

B. Comparison with the variational approximation 

Comparison of the predictions of the VA with numerical results in two typical cases 
(for stable and decaying solitons, corresponding to = 5.4 and = 4.6, respectively) is 
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FIG. 4: (Color online) Comparison of the evolution of the soliton's amplitude, as predicted by 
the variational approximation ( "VA" ) and found from direct simulations of the partial differential 
equation ^ ("PDE"), for N = 5.4 (a stable soliton) and N = 4.6 (a decaying state). 

presented in Fig. |4j The discrepancy in the dependence of the soliton's amplitude on time, 
observed in the former case, is a known feature [T71 |38l HI], explained by the fact that the 
oscillations predicted by the VA are damped by the radiation loss, which is not taken into 
regard by the VA based on simple ansatz (|3]). A modification of the VA aimed to account 
for the emission of radiation, in some approximation, was proposed but it is quite 
cumbersome even in the ID setting. 

The VA predicts eigenfrequency uq of intrinsic oscillations of the weakly perturbed 2D 
soliton trapped in the static OL, see Eq. Because this feature may be important to 

understand the response of the soliton to the periodic time modulation of the lattice, in 
terms of possible resonances (see below), it is necessary to check whether the presence of the 
eigenfrequency is confirmed by simulations of full equation (IT]) with the static lattice, i.e.. 
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FIG. 5: (Color online) The power spectrum of small random perturbations around the stable 
soliton trapped in the static lattice (e = 0), for Vq = 0.25 and N = 5.31. 

e = 0. To this end, in Fig. [5] we display the power spectrum of small oscillations around 
the soliton caused by the addition of a small random perturbation to it, for Vq = 0.25 
and = 5.31, which corresponds to 1 — A^/A^max ~ 0.155. As mentioned above, in this 
case Eqs. ([s]) and ^ predict the eigenfrequency to be cuq ^ 1.35. The spectrum in Fig. 
|5] clearly shows the main peak quite close to this point. Peaks corresponding to higher- 
order resonances (see above) can also be recognized in the figure. To produce the spectrum 
shown in FigjS} we simulated the evolution of the soliton up to a very long time, t = 1200, 
eliminating a contribution from a relatively short initial stage, which featured a transient 
behavior. 
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C. Stability diagrams 



Results produced by the systematic numerical analysis of the stability of the 2D soliton 
under the action of the "lattice management" are collected in Fig. [6|^a), which displays the 
stability region in the plane of the modulation parameters, uj and e, for fixed Vq = 0.25 and 

= 5.905. Note that this norm again (like in the cases of Figs. [T]and[2]) slightly exceeds 
the collapse threshold in the static model, which is given by Eq. ([T]). The stability region 
in the plane of e and A^, for the same OL strength, Vq = 0.25, as in Fig. |6j and = 4, is 
displayed in Fig. [6|^b). The latter plot explicitly demonstrates the growth of the collapse 
threshold, A^maxj with the increase of the modulation amplitude. Note that not the entire 
parameter region below the stability border in Fig. [6]^b) corresponds to stable solitons; if N 
is too small, the solitons are unstable against decay, see below. 

In Fig. |6](a), the horizontal line denotes the maximum possible value of the modulation 
amplitude in the model of the time-modulated OL, £max = 2. A noteworthy fact is that the 
stability region can extend up to this limit, i.e., 100% modulation depth (recall that, in the 
case when the stabilization of 2D solitons was provided by the quasi-lD lattice, the stability 
limit corresponded to shallow modulation [3^). The results for e > 2 are included too in 
Fig. [6]^a), as they may find a different physical realization, corresponding to a superposition 
of moving OLs [33J (as mentioned above). 

A feature obvious in Fig. [6]^a) is a resonant-like dependence of the stability border on 
UJ. This may be a manifestation of the fundamental resonance at uj close to ujq and higher- 
order resonances at multiple values of uj. As argued above, the resonance may help to arrest 
the collapse, by forcing the vibrating soliton to spend less time in the state where it is 
"dangerously" narrow. Since the value of the norm corresponding to Fig. |6](a) is close to 
A^max, one may expect that the corresponding fundamental-resonance frequency should be 



close to one given by Eq. (10), i.e., ujq = 2, for Vq = 0.25 [the same is given by Eqs. 



and ([9j) in the limit of N/ A^max 0] . Indeed, the picture in Fig. |6j^a) is consistent with the 
expectation of the fundamental and higher-order resonances aX uj = nujo, n = 1,2,3,4 (the 
picture also suggests the existence of a very strong resonance close to lj = 5ujq, but that one 
falls deeply into the unphysical region, e > 2). 

As mentioned above, the soliton whose norm is too small cannot stabilize itself and decays 
into quasi-linear waves. Therefore, along with the upper stability boundary, A^ = A^max, that 
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FIG. 6: (Color online) (a) The stability region in the plane of the modulation parameters, 00 
and e, for Vq = 0.25 and N = 5.905. The area relevant to the model with the periodically 
modulated optical lattice corresponds to e < 2, where the time-dependent amplitude in Eq. ([T]), 
1 + (e/2) cos(Ljt), does not change its sign, (b) The collapse threshold versus the modulation 
amplitude, e, for Vq = 0.25 and uj = A. 

guarantees the absence of the collapse, it is necessary to identify the lower one, = A/'min, 
which secures the stability of the soliton against the decay. Both boundaries, as predicted 
by the VA and found from the direct simulations, are displayed in Fig. [7], in the form of 
dependences A^'maxlK)) and A^min(Vo), at fixed e = 0.5 for several different values of u. In 
each panel, the stability region is A^min < N < A^max- 

We observe in Fig. [Tjthat dependences A^min(Vo) produced by the VA are in reasonable 
agreement with the results of direct simulations of the GPE for modulation frequencies 
uj > 2. However, there is a conspicuous discrepancy between the VA and GPE for uj close to 
1. On the other hand, this range features strong resonances in the perturbation spectrum, see 
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FIG. 7: (Color online) The dependence of the stability boundaries, A^max and A^^mim on the OL 
strength, Vq, at e = 0.5 and different values of the modulation frequency, lo. Solid and dashed lines 
show, severally, A'max and A^min as found from the simulations, while the dotted and dashed-dotted 
lines represent, respectively, A'^max and N^i^ as predicted by the variational approximation. The 
arrows in (c) and (d) indicate values A^'max for the stationary case, e = 0. 

Fig. |5} Therefore, it is interesting to explore the A^min(Vo) dependence in region 0.5 < to < 2, 
where the resonances may lead to decay of the soliton. These results are displayed in Fig. 
[sj which shows that plots A^min(Vo) essentially differ, in this region, from their counterparts 
both in the static model (see the curve for e = in Fig. |8]) and in the high-frequency region, 
uj > 2: actually, A^min increases due to the resonant decay of the solitons. The fact that the 
resonances are narrow enough, as seen in Fig. |5| may explain a non-monotonous character 
of dependences A^min(Vo) for u = 1.10 and 1.25. 

Further, dependences of A^max and A^min on the modulation frequency are displayed in 
Fig. |9| for the same modulation amplitude as above, e = 0.5, and several different values of 
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FIG. 8: (Color online) Dependencies N^i^iVo) for e = 0.5 and different modulation frequencies co 
in the resonant area, 0.5 < u < 2 (see details in text). The curve pertaining to e = (the static 
model) is included for comparison. 

Vq. It is quite natural that dependences N^i^{u), as predicted by the VA and generated by 
direct simulations of Eq. ([T]), are close to each other for smaller values of the OL strength, 
Vo = 0.1 and 0.25, while at Vq = 0.5 the discrepancy between them is considerable. Note 
non-monotonous behavior of A^min(ti^) in resonant zone 0.5 < u < 2, which correlates to 
peculiarities of dependences A^min(Vo) in the same zone, cf. Fig. |8] 

V. CONCLUSIONS 

We have studied the dynamics of 2D solitons in the model of BEC trapped in the square- 
shaped OL (optical lattice) whose strength it subject to the periodic time modulation. Being 
quite feasible for the experimental implementation, the model belongs to a broad class of 
schemes of the periodic management of solitons [23] . By means of the VA (variational ap- 
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FIG. 9: (Color online) The dependence of the upper (a) and lower (b) stability boundaries, A^'max 
and N^in, on modulation frequency uj for e = 0.5 and different values of OL strength Vq, which 
are specified in the box. Labels "v" and "g" pertain, severally, to the curves predicted by the 
variational approximation and those found from direct simulations of Gross-Pitaevskii equation 
(the VA predicts flat value iVmax = 27r, see text). The arrows identify those curves which may 
seem indisnguishable in the black-and-white rendition of the figure. 

proximation) and direct systematic simulations of the underlying Gross-Pitaevskii equation, 
we have identified stability regions for the solitons in the parameter space of the model, 
including both maximum and minimum values of the norm, as summarized in Figs. [6} [7], [8] 
and |9} A remarkable feature demonstrated by these results is that the stability limit may 
reach the maximum (100%) modulation depth. It is noteworthy too that an increase of the 
collapse threshold in predicted in comparison with its classical value in the static situation, 
which corresponds to the norm of the Townes soliton. The stability borders predicted by 
the VA are found to be in reasonable agreement with the numerical results. In the plane of 



18 



the modulation frequency and amplitude, the stability boundary features a salient resonant 
structure, which may also be qualitatively explained by means of the VA. In particular, the 
minimum norm necessary for the stability of the soliton demonstrates strong variations in 
the region of low modulation frequencies, due to the possibility of the decay of the soliton 
enhanced by the resonance. 

The analysis reported in this work can be extended in several directions. It may be 
interesting to study interactions between the solitons in this setting, and identify stability 
limits for vortex solitons, as well as for 2D gap solitons in the model combining the repulsive 
nonlinearity and lattice management. Another extension may be made in the direction of 
an anisotropic lattice management, i.e., applying time modulations shifted by n to the two 
ID sublattices. 
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